#ifndef UTILS_CC
#define UTILS_CC

#include<string>
#include<fstream>
#include<sstream>
#include<iostream>
#include<algorithm>
#include<iterator>
#include<vector>
#include<map>
#include<cstdlib>
#include<ctime>
#include<cmath>
#include<random> // this is for C++ 0x
#include "def.hpp"
#include "utils.hpp"
#include "lngamma.hpp"
#include "digamma.hpp"

using namespace std;

namespace Puza{
  double Utils::get_sum(double*& buffer, int& K){
    double result = 0.0;
    for(int k=0; k < K; k++){
      result += buffer[k];
    }
    return result;
  }
  
  double Utils::log_sum(double log_a, double log_b){
    if(log_a < log_b)
      return log_b + log(1 + exp(log_a - log_b));
    else
      return log_a + log(1 + exp(log_b - log_a));
  }

  double Utils::log_sub(double log_a, double log_b){
    if(log_a < log_b)
      return log_b + log(1 - exp(log_a - log_b));
    else
      return log_a + log(1 - exp(log_b - log_a));
  }
  
  double Utils::log_sum_array(double*& buffer, int K){
    // first find the maximum
    double max_x = buffer[0];
    double result = 0.0;  
    for(int k=0; k < K; k++){
      if(max_x < buffer[k])
	max_x = buffer[k];
    }
    // update
    for(int k=0; k < K; k++){
      result = result + exp(buffer[k] - max_x);
    }
    result = log(result) + max_x;
    return result;
  }
  
  double Utils::log_sum_vector(vector<double>& buffer){
    // first find the maximum
    double max_x = buffer[0];
    double result = 0.0;
    for(unsigned int k=0; k < buffer.size(); k++){
      if(max_x < buffer[k])
	max_x = buffer[k];
    }
    // update
    for(unsigned int k=0; k < buffer.size(); k++){
      result = result + exp(buffer[k] - max_x);
    }
    result = log(result) + max_x;
    return result;
  }
  
  
  double Utils::digamma(double x){
    int return_value;
    return PuzaThirdParty::digama(x,&return_value);
  }
  
  double Utils::log_gamma(double x){
    int return_value;
    return PuzaThirdParty::lngamma(x,&return_value);
  }
  
  double Utils::sum_log_gamma(double*& buffer, int& K){
    double s = 0.0;
    for(int i=0; i < K; i++){
      s = s + log_gamma(buffer[i]);
    }
    return s;
  }
  
  double Utils::l1_distance(double*& a, double*& b, int& K){
    double result = 0.0;
    for(int k=0; k < K; k ++){
      result = result + fabs(a[k] - b[k]);
    }
    return result;
  }
  
  double Utils::l2_distance(double*& a, double*& b, int& K){
    double result = 0.0;
    for(int k=0; k < K ; k++){
      result = result + (a[k] - b[k]) * (a[k] - b[k]);
    }
    result = sqrt(result);
    return result;
  }
  
  double Utils::Min(double a, double b){
    if(a < b)
      return a;
    else
      return b;
  }
  
  double Utils::Max(double a, double b){
    if(a > b)
      return a;
    else
      return b;
  }
  
  double Utils::Trunc(double vj, double alpha){
    if(vj <=0)
      return Min(0,vj + alpha);
    else
      return Max(0,vj - alpha);
  }
  
  int Utils::uniform_int_sample(int lower,int upper){
    uniform_int_distribution<int> uniform_int(lower,upper);
    return uniform_int(PuzaDEF::Instance()->random_eng);
  }
  
  double Utils::uniform_double_sample(double lower,double upper){
    uniform_real_distribution<double> uniform_real(lower,upper);
    return uniform_real(PuzaDEF::Instance()->random_eng);
  }
};
#endif
